% Optimal Harp Seal Harvest with 2 Stochastic Variables
% Code adapted from demsc02 in CompEcon Toolbox by Charles Sims, University of Tennessee
% adapted July 2019

% This code (along with the files listed below) solve the social planner's
% optimal harp seal harvest.  The harp seal population (S) is stochastic and
% dependent on sea ice level.  Mean sea ice (Z) is declining over time
% the 
% dS=[rS(1-S/K)-q]dt+sig(S)*dw
% dZ=-gamma*Z*dt

% SOURCES:
% For more information on the value function approximation procedure see: 
% Judd, K. 1998. Numerical Methods in Economics. Cambridge, MA: MIT Press.
% Miranda, M., and P.L. Fackler. 2002. Applied Computational Economics and Finance: The MIT Press

%For more inforamtion on the stochastic harvest model see:
% Sims, C., R. Horan, B. Meadows. 2018. Come on feel the noise: ecological foundations in stochastic bioeconomic models. Natural Resource Modeling
% Horan and Sims. Managing Harp Seals After the Cold November Rain: Accounting for Arctic Warming and Ecological Stochasticity. Working paper.

% DEPENDENCIES:
% The following files are needed to run this code.  Make sure they are added to Matlab's "path" so that Matlab can find these files
% mff2varharp12.m (model function file: provides model structure and functional forms)
% CompEcon Toolbox
% scsolve_sigx2.m (variant of scsolve from CompEcon that allows for the control variable in the diffusion term and calculate the 2nd derivative of value function)
% itsimul2.m (variant of itosimul.m from CompEcon that allows model function file with Vss)

% NOTES:
% This code uses two scaling procedures.  The first is to scale S so that
% it is of the same relative magnitude as Z (called scale below).  The second is to scale the
% benefits function so it is closer to 1 (called scale2 below).  See scaling the seal model
% 2-17.docx for more information.

%close all
optset('scsolve','defaults');
disp(' ')
disp('Renewable Resource Management with Stochastic Growth and Climate Change')

%% Parameter values
% SET PARAMETER SCENARIO
% NOTES: Results for scenario 0 and scenario 2 are included in the JAERE
% paper. Scenario 5 is inclued in the appendix. Run scenario 1 first to
% create Figure 5 from JAERE paper.
% scenario=0 Baseline (type III; mean variance habitat loss): climate change that increases sea ice variability
% scenario=1 No climate change: sea ice is not declining (gamma=0)
% scenario=2 climate change doesn't affect sea ice volatility (type I; mean-only habitat loss): sig_z=sig0-(theta*sig1*Z+(1-theta)*sig1*Zmax)
% scenario=3 no environmental variance (gamma>0; sig0=sig1=0)
% scenario=4 no demographic variance (D=0)
% scenario=5 climate change only increases variability (type II; variance-only habitat loss) (gamma=0 and sig0 increases)
percent=1; % how to report differences: percent or level values (percent=1 for percent)
scenario=0; % choose scenario above
benchmark=1; % choose which scenario is your benchmark when calculating Figure 5 from JAERE paper

% ecological parameters (S=harp seal stock) 
% modified Schoener (1973) logistic growth: dS/dt=S[rho(1-lambda*S)-zeta]=rS(1-S/K) where r=(rho-zeta) and K=(rho-zeta)/(rho*lambda)
nu      = 0.0016; % bycatch rate
zeta  = 0.0253+nu; % adult natural mortality rate
lambda = 7.74E-8; % density dependent cost that reduces recruitment
scale     = 1000000; % scalar for harp seal population; scale=1000000 means S is in millions of harp seals
K1    = (1/lambda)/scale;
b     = 0.165; % birth rate 
m0    = 0.86; % intercept of juvenile mortality prior to density dependent impacts
m1    = 0.073; % slope of juvenile mortality prior to density dependent impacts

% economic parameters
% linear demand fxn for harvest: p(qS)=p0-p1*q*S
% marginal cost of harvest if function of current sea ice, z: c(z)=c0-c1*z
p0=20.99; % intercept of demand function
p1=3.39E-5; % slope of demand function
Qmax=p0/p1; % Price<0 if total harvest exceeds Qmax
c0=15.407*1.25; %intercept of marginal cost function
c1=0.673*1.25; % slope of marginal cost function
u0=126026.85; % parameter from existence value function
delta   = 0.05; % discount rate  
scale2  = 1000000^2; % scale for harp seal soclal net benefits; scale2=1000000 means B and V are in millions of $

% Benchmark sea ice parameters (Z=mean sea ice)
sig0  =38.2; % intercept of var{z|Z}
sig1  =3.24; % slope of var{z|Z}
gamma =0.019; % rate of sea ice decline 

if scenario==1 % no climate change
    gamma =0; % rate of sea ice decline
end
if scenario==2 % climate change doesn't increase variability
    theta = 0; % parameter that determines whether sigma_z is a function of Z
else
    theta = 1; % parameter that determines whether sigma_z is a function of Z
end
if scenario==3 % determinisitc climate change
    sig0  =0; % intercept of var{z|Z}
    sig1  =0; % slope of var{z|Z}
    gamma=0;
end
if scenario==4 % dummy variable turns dem variance on/off (=1 on; =0 off)
    D=0;
else
    D=1;
end
if scenario==5
    gamma = 0;  % rate of sea ice decline
    sig0 = 2*38.2; % intercept of var{z|Z}
end

Zdet  =sig0/sig1; % variance of z and environmental variance is 0 at Zdet (Z>Zdet implies negative variance; Z=Zdet -> demographic variance only)

% Impose boundary condition on approx value function for generating green diff plots in paper
bound=0; %if bound=1, value function at Z=0 in all scenarios will be the same as no climate change scenario 1

font=16;

%% Set up the approximation space and define the number of approximation nodes
% Solves model over nondimensional state space.  
Smin=1; Smax=K1;           % stock limits in millions of seals
Zmin=0; Zmax=8;   % mean sea ice limits
Snodes=50; Znodes=50; 
x1=linspace(Smin,Smax,Snodes*10+1);  % 10 is the measure for nres in scsolve
Sstep=x1(2)-x1(1); % step size along S dimension in millions
x2=linspace(Zmin,Zmax,Znodes*10+1);
Zstep=x2(2)-x2(1);

%% Set model structure
clear model
model.func='mff2varharp_final';
model.params={p0,p1,c0,c1,u0,delta,lambda,zeta,K1,b,m0,m1,gamma,sig0,sig1,scale,scale2,theta,D,Zmax};
model.a=[0 0];  % lower bounds for the stochastic processes S and Z  
model.b=[inf inf]; % upper bounds for the stochastic processes S and Z  

%% Choose upwind or centered finite difference derivatives
% Upwinding method and piecewise linear are highly recommended for convergence issues.
upwind=1;                   % 1 for upwind, 2 for centered finite difference
family=1;                   % 1 for piecewise linear, 2 for cubic spline

if (family==1),             % define approximaion family 
    af='lin'; 
else 
    af='spli'; 
end

fspace=fundefn(af,[Snodes Znodes],[Smin Zmin],[Smax Zmax]);

    snodes=funnode(fspace);     % define the nodes in state space (dimension=2)
    s=gridmake(snodes);         % create the grid with the nodal points
                                % S=S(:,1), Z=S(:,2)   

                                
%% Call the stochastic control solver
xinit=0.05; % initial guess for q; data suggests harvests have been around 5% of the population 
% scsolve_sigx2 adjusts scsolve to allow for the harvest rate in the
% variance term of the SDE and for Vss to be an input in the calculation of the harvest rate 
[cv,S,v,x,resid] = scsolve_sigx2(model,fspace,snodes,[],xinit);

% Impose boundary condition for V at Z=0
if bound==1
cvmat=reshape(cv,[Snodes,Znodes]);
if scenario==benchmark
    cvbench=cv(1:Snodes);
end
if scenario~=benchmark
    cv=cv(Snodes+1:end); % remove first Snodes elements
    cv=[cvbench; cv]; % add coefficients from benchmark when Z=0
end
end

%Evaluate the value fucntion with boundary condition imposed
v=funeval(cv,fspace,S,[0,0])*scale2; % rescale the value function (now in terms of seals)
v=reshape(v,[size(x1,2),size(x2,2)]); % convert value function vector to matrix 

% Calculate partial derivatives of value function wrt S
% NOTE: the calculation of derivatives uses a linearization procedure that 
% means the scaling of these derivatives differs from what is described in scaling the seal model 2-17.docx  
Vsa=funeval(cv,fspace,S,[1,0])*scale2/scale; % evaluate first derivative of approximated value function
Vsa=reshape(Vsa,[size(x1,2),size(x2,2)]); % approx marginal shadow value of harp seals ($/seal)
Vssa=funeval(cv,fspace,S,[2,0])*scale2/scale; % evaluate second derivative of approximated value function
Vssa=reshape(Vssa,[size(x1,2),size(x2,2)]);
Vsssa=funeval(cv,fspace,S,[3,0])*scale2/scale; % evaluate third derivative of approximated value function
Vsssa=reshape(Vsssa,[size(x1,2),size(x2,2)]);
% Calculate partial derivatives of value function wrt Z
Vza=funeval(cv,fspace,S,[0,1])*scale2; % evaluate first derivative of approximated value function
Vza=reshape(Vza,[size(x1,2),size(x2,2)]);
Vzza=funeval(cv,fspace,S,[0,2])*scale2; % evaluate second derivative of approximated value function
Vzza=reshape(Vzza,[size(x1,2),size(x2,2)]);
% Calculate cross partials
Vzsa=funeval(cv,fspace,S,[1,1])*scale2/scale; % evaluate first derivative of approximated value function
Vzsa=reshape(Vzsa,[size(x1,2),size(x2,2)]);
Vzssa=funeval(cv,fspace,S,[2,1])*scale2/scale; % evaluate first derivative of approximated value function
Vzssa=reshape(Vzssa,[size(x1,2),size(x2,2)]);

% Calculate optimal harvest and other metrics
[Smat, Zmat]=ndgrid(x1,x2);  
Qmat=x.*Smat.*scale; % total harvest
Qmat_scaled=x.*Smat; % total harvest in millions of seals
Price=p0-p1*Qmat; % $ per harp seal
Cost=c0-c1*Zmat; % marginal cost function
r=b*(1-(m0-m1*Zmat))-zeta;
K=(1-zeta./(b*(1-m0+m1*Zmat)))*K1;
K_Z=zeta*m1./(b*(1-m0+m1*Zmat).^2)*K1; % derivative of K wrt Z
Kdummy=Smat>K;
Growth=r.*Smat.*scale.*(1-Smat./K);
marg_Growth=r.*(1-2*Smat./K);
NetGrowth=Growth-Qmat;
env_var=(b^2*m1^2*(sig0-(theta*sig1*Zmat+(1-theta)*sig1*Zmax)).*(1-Smat/K1).^2).*(Smat*scale).^2; % environmental variance in seals
dem_var=D*(b*(1-m0+m1*Zmat).*(1-Smat/K1)+zeta+x).*Smat*scale; % demographic variance in seals
var_dif=env_var-dem_var;
Sigma2=env_var+dem_var; % variance in seals
Sigma=Sigma2.^0.5; % sigma term in seals
Sigma_scaled=(env_var/scale^2+dem_var/scale).^0.5; % sigma term in millions of seals
Existence=u0*log(Smat*scale); % harp seal existence value
Profit=Qmat.*p0-0.5*p1.*(Qmat).^2-Cost.*Qmat; % Profit
Rent=p0-p1*Qmat-Cost; % marginal resource rent
SSFmat=(Profit+Existence); % social surplus function: int(p(w)dw)-c(z)q*S+u(s)

% Calculate risk terms for S adjoint
% do everything with scaled S (millions of S) for consistency -> this means multiply Vs by scale below
As=-Vssa./Vsa; % index of absolute risk aversion (don't need to scale this since divided by two partials)
RPs=As.*(0.5*(1-2*Smat/K1).*(D*b*(1-m0+m1*Zmat)+b^2*m1^2*(sig0-(theta*sig1*Zmat+(1-theta)*sig1*Zmax)).*2*Smat.*(1-Smat/K1))+D*zeta); % LHS risk premium
CGRAs=0.5.*Sigma_scaled.^2.*(Vsssa./Vsa); % RHS risk adjustment term
DetCG=(Vssa./Vsa).*(r.*Smat.*(1-Smat./K)-x.*Smat)-Vzsa./(Vsa*scale)*gamma.*Zmat; % Deterministic captial gains term
Ret=(marg_Growth)+(u0./(Smat*scale))./(Vsa*scale)+DetCG+CGRAs; % expected return to conservation
CG_Ratio=CGRAs./Ret; % Percent of expected return due to risk adjustment
NetSrisk=RPs-CGRAs;
adjointS_stoch=delta+RPs-Ret; % adjoint equation for S with risk corrections
adjointS_det=delta-(marg_Growth)-DetCG; % adjoint equation for S without risk corrections

% Calculate risk terms for Z adjoint (do everything with scaled S for consistency)
Az=-Vzsa.*scale./Vza; % index of absolute risk aversion (scale numerator since a partial wrt S)
RPz=Az*0.5*b*m1.*Smat.*(1-Smat./K1).*(D*1-b*m1*sig1*Smat.*(1-Smat./K1)); % LHS risk premium
CGRAz=0.5*Sigma_scaled.^2.*(Vzssa.*scale./Vza); % capital gains risk adjustment for the Z adjoint eqn
Loss=-(-c1*x.*Smat)./Vza+(b*m1*Smat.*(1-Smat./K)+r.*(Smat./K).^2.*K_Z).*(Vsa.*scale./Vza)-gamma+(r.*Smat.*(1-Smat./K)-x.*Smat).*(Vzsa.*scale./Vza)-(Vzza./Vza).*gamma.*Zmat+CGRAz;
omega=Loss-RPz;
adjointZ_stoch=delta+RPz-Loss; % adjoint equation for Z with risk corrections
comb_adj=omega+RPs-Ret; % combined adjoint condiiton

% harvest equation adjustment term (May 12, 2021)
Ret_adj=(marg_Growth)+(u0./(Smat*scale))./(Vsa.*scale)+CGRAs; % expected return to conservation (note DetCG disappears since this is evaluated for E{dS}=Zdot=0)
Loss_adj=-(-c1*x.*Smat)./Vza+(b*m1*Smat.*(1-Smat./K)+r.*(Smat./K).^2.*K_Z).*(Vsa.*scale./Vza)-gamma-(Vzza./Vza).*gamma.*Zmat+CGRAz; % note that several terms are missing relative to the loss equation above because this is evaluated for E{dS}=Zdot=0 
adjterm=Vsa./(-2*Vssa).*(delta+RPs-Ret_adj+((delta+RPs-Ret_adj+(-c1)*gamma*Zmat./(Vsa.*scale)).^2+4*(Vza*gamma.*Zmat./(Vsa.*scale))*(-Vssa./Vsa)*(Loss_adj-delta-RPz)).^0.5);
adjterm_nochange=Vsa./(-Vssa).*(delta+RPs-Ret_adj);
adjterm_approx=Qmat_scaled-(r.*Smat.*(1-Smat./K));

% fit surface for policy function
Svec=Smat(:);
Zvec=Zmat(:);
xvec=x(:);
xfun=fit([Svec,Zvec],xvec,'poly31');
xfit=feval(xfun,Smat,Zmat);
fiterror=(x-xfit);
fitMSE=sum(sum(fiterror.^2))/numel(xfit);

% Calculate changes in value function and harvest rate under differnt scenarios
if scenario==benchmark
vbench=v; % set benchmark value fucntion
qbench=x; %set benchmark harvest rate
Qmatbench=Qmat; % set benchmark total harvest
SRentbench=Vsa; % set benchmark rent
ZRentbench=Vza; % set benchmark rent
adjbench=-NetGrowth; % set benchmark adjustment term in expression 15
end

if scenario~=benchmark
    if percent==1
        vadj=v-vbench;
        vadj=max(vadj(1,:)); % shifts benchmark value function.  Used to avoid increase in value function from habitat loss.  Need to change depending on Type.
        vchange=(v-(vbench+vadj))./(vbench+vadj); % calculate percent change in value function is considering alternative
        qchange=(x-qbench)./(qbench+1E-14); % calculate percent change in harvest rate is considering alternative
        Qchange=(Qmat-Qmatbench)./(Qmatbench+1E-14); % calculate percent change in harvest is considering alternative
        SRentchange=(Vsa-SRentbench)./SRentbench; 
        ZRentchange=(Vza-ZRentbench)./ZRentbench; 
        
    else
        vchange=v-vbench; % calculate level change in value function
        qchange=x-qbench; % calculate level change in harvest rate
        Qchange=Qmat-Qmatbench; % calculate level change in harvest
        SRentchange=Vsa-SRentbench; 
        ZRentchange=Vza-ZRentbench;
        
    end
end

%% Figure 3 Plots From JAERE Paper
figure(1) % growth rate and carrying capacity
yyaxis left
plot(Zmat(1,:),r(1,:));
ylabel('intrinsic growth rate, r','fontsize',18);
yyaxis right
plot(Zmat(1,:),K(1,:)*1000000);
ylabel('carrying capacity, K','fontsize',18);
xlabel('mean sea ice (Z)','fontsize',18);

figure(2) % growth function
hold on
%mesh(Smat*scale,Zmat,Growth);
[C,h]=contour(Smat*scale,Zmat,Growth);
clabel(C,h,'FontSize',14)
contour(Smat*scale,Zmat,Kdummy,[1,1],'LineWidth',2,'LineColor',[0 0 0],'LineStyle','--');
xlabel('harp seals (S)','FontSize',18)
ylabel('mean sea ice (Z)','FontSize',18)
%zlabel('Net natural seal growth, G(S,\rho)')
title('Mean net natural seal growth, G(S,Z)','FontSize',18)

figure(3) % plot population variance
hold on
%mesh(Smat*scale,Zmat,Sigma2);
contourf(Smat*scale,Zmat,RPs,[0,0],'LineColor','none')
%[C,h]=contour(Smat*scale,Zmat,Sigma2,[1e10 2e10 3e10 4e10 5e10],'LineColor',[0 0 0]);
%clabel(C,h,'FontSize',14);
contour(Smat*scale,Zmat,Sigma2,[1e10 2e10 3e10 4e10 5e10],'LineColor',[0 0 0])
%contour(Smat*scale,Zmat,RPs,[0,0],'ShowText','on','LineWidth',2,'LineColor',[0 0 0])
contour(Smat*scale,Zmat,Kdummy,[1,1],'LineWidth',2,'LineColor',[0 0 0],'LineStyle','--')
colormap(gray)
caxis([-0.5,0.5])
xlabel('harp seals (S)','FontSize',18)
ylabel('mean sea ice (Z)','FontSize',18)
%zlabel('variance of seal population changes, \sigma(S,Z,q)^2')

figure(4) % plot demographic variance
hold on
%mesh(Smat*scale,Zmat,dem_var);
%[C,h]=contour(Smat*scale,Zmat,dem_var);
%clabel(C,h,'FontSize',14)
contour(Smat*scale,Zmat,dem_var)
contour(Smat*scale,Zmat,Kdummy,[1,1],'LineWidth',2,'LineColor',[0 0 0],'LineStyle','--')
xlabel('harp seals (S)','FontSize',18)
ylabel('mean sea ice (Z)','FontSize',18)
%zlabel('demographic variance')

figure(5) % plot environmental variance
hold on
%mesh(Smat*scale,Zmat,env_var);
%[C,h]=contour(Smat*scale,Zmat,env_var,[1e10 2e10 3e10 4e10 5e10]);
%clabel(C,h,'FontSize',14);
contour(Smat*scale,Zmat,env_var,[1e10 2e10 3e10 4e10 5e10])
contour(Smat*scale,Zmat,Kdummy,[1,1],'LineWidth',2,'LineColor',[0 0 0],'LineStyle','--')
xlabel('harp seals (S)','FontSize',18)
ylabel('mean sea ice (Z)','FontSize',18)
%zlabel('environmental variance')

%% Figure 4 Plots From JAERE Paper

figure(6);  % value function
hold on
%mesh(Smat*scale,Zmat,v/1000000)
[C,h]=contour(Smat*scale,Zmat,v/1000000);
clabel(C,h,'FontSize',14)
contour(Smat*scale,Zmat,Kdummy,[1,1],'LineWidth',2,'LineColor',[0 0 0],'LineStyle','--');
xlabel('harp seals (S)','FontSize',18)
ylabel('mean sea ice (Z)','FontSize',18)
%zlabel('Value Function, V (million $)')
title('Value Function, V (million $)','FontSize',18)

figure(7); % plot total harvest (q*S) 
hold on
%mesh(Smat*scale,Zmat,Qmat)
contourf(Smat*scale,Zmat,-Qmat+Growth,[0,0],'LineColor','none')
[C,h]=contour(Smat*scale,Zmat,Qmat,'LineColor',[0 0 0]);
clabel(C,h,'FontSize',14)
contour(Smat*scale,Zmat,Kdummy,[1,1],'LineWidth',2,'LineColor',[0 0 0],'LineStyle','--');
colormap(gray)
caxis([-0.5,0.5])
xlabel('harp seals (S)','FontSize',18)
ylabel('mean sea ice (Z)','FontSize',18)
%zlabel('total harvest (q^*S)')
title('total harvest, q^*(S,Z)S','FontSize',18)

figure(8); % plot shadow value of harp seals 
hold on;
%mesh(Smat*scale,Zmat,Vsa)
[C,h]=contour(Smat*scale,Zmat,Vsa);
clabel(C,h,'FontSize',14)
contour(Smat*scale,Zmat,Kdummy,[1,1],'LineWidth',2,'LineColor',[0 0 0],'LineStyle','--')
xlabel('harp seals (S)','FontSize',18)
ylabel('mean sea ice (Z)','FontSize',18)
%zlabel('$/seal')
title('$/seal','FontSize',18)

figure(9); % plot shadow value of sea ice 
hold on;
%mesh(Smat*scale,Zmat,Vza/1000000)
[C,h]=contour(Smat*scale,Zmat,Vza/1000000);
clabel(C,h,'FontSize',14)
contour(Smat*scale,Zmat,Kdummy,[1,1],'LineWidth',2,'LineColor',[0 0 0],'LineStyle','--')
xlabel('harp seals (S)','FontSize',18)
ylabel('mean sea ice (Z)','FontSize',18)
%zlabel('million $/percent area of ice cover') 
title('million $ per percent area of ice cover','FontSize',18)

%% Figure 5 Plots from JAERE Paper
if scenario~=benchmark
figure(10) % plot change in the value of fishery
hold on
contourf(Smat*scale,Zmat,vchange,[0,0],'LineColor','none')
contour(Smat*scale,Zmat,vchange,'ShowText','on','LineColor',[0 0 0])
contour(Smat*scale,Zmat,vchange,[0,0],'ShowText','on','LineWidth',2,'LineColor',[0 0 0])
contour(Smat*scale,Zmat,Kdummy,[1,1],'LineWidth',2,'LineColor',[0 0 0],'LineStyle','--')
colormap(summer)
caxis([-0.5,0.5])
xlabel('harp seals (S)')
ylabel('mean sea ice (Z)')
%title('change in value of fishery')

figure(11) % plot change in the harvest
hold on
contourf(Smat*scale,Zmat,Qchange,[0,0],'LineColor','none')
contour(Smat*scale,Zmat,Qchange,'ShowText','on','LineColor',[0 0 0])
contour(Smat*scale,Zmat,Qchange,[0,0],'ShowText','on','LineWidth',2,'LineColor',[0 0 0])
contour(Smat*scale,Zmat,Kdummy,[1,1],'LineWidth',2,'LineColor',[0 0 0],'LineStyle','--')
colormap(summer)
caxis([-0.5,0.5])
xlabel('harp seals (S)')
ylabel('mean sea ice (Z)')
%title('change in harvest')
end





  